o 

(N 

in 



> 



This manuscript is accepted by RAA 

Non-similar collapse of singular isothermal spherical molecular 
cloud cores with nonzero initial velocities 

Mohsen Nejad-Asghar^'^'^ 

^Department of Physics, University of Mazandaran, Baholsar, Iran 

"^Research Institute for Astronomy and Astrophysics of Maragha, Maragha, Iran 

^Department of Physics, Damghan University, Damghan, Iran 



< 

^ , ne j adasgharSumz . ac . ir 



ABSTRACT 



6 

c^ I Theoretically, stars have been formed from the collapse of cores in the molec- 



ular clouds. Historically, the core had been assumed as an singular isothermal 

sphere (SIS), and the collapse had been investigated by a self-similar manner. 

lO ! This is while the rotation and magnetic field lead to non-symmetric collapse so 

^ I that a spheroid shape may be occurred. Here, the resultant of the centrifugal 

O I force and magnetic field gradient is assumed to be in the normal direction of the 

^ I rotational axis, and its components are supposed to be a fraction /3 of the local 

O I gravitational force. In this research, a collapsing SIS core is considered to find the 

importance of the parameter /3 for oblateness of the mass shells which are above 

the head of the expansion wave. We apply the Adomian decomposition method 

r> ' to solve the system of nonlinear partial differential equations because the collapse 

c^ I does not occur in a spherical symmetry with self-similar behavior. In this way, 

we obtain a semi-analytical relation for the mass infall rate M of the shells at 

the envelope. Near the rotational axis, the M decreases with increasing of the 

non-dimensional radius C,, while a direct relation is observed between M and ^ 

in the equatorial regions. Also, the values of M in the polar regions are greater 

than the equatorial values, and this difference is more often at smaller values of 

^. Overall, the results show that before reaching the head of expansion wave, 

the visible shape of the molecular cloud cores can evolve to oblate spheroids. 

The ratio of major to minor axes of oblate cores increases with increasing the 

parameter (3, and its value can approach to the apparently observed elongated 

shapes of cores in the maps of molecular clouds such as Taurus and Perseus. 
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Subject headings: ISM: clouds - ISM: evolution - star: formation - methods: 
analytical 



1. Introduction 

Great deals are now known about dense cores in the molecular clouds that are the 
progenitors of protostars (e.g., di Francesco et al. 2007, Ward-Thompson et al. 2007). We 
know that, approximately, all cores in the maps of molecular clouds seem apparently to be 
elongated rather than spherical. For example, we can refer to the recent work of Curtis and 
Richer (2010) for two-dimensional ellipticity of cores in the Perseus molecular cloud, or to 
the old report of Myers et al. (1991) for apparent elongated shapes of cores in the Taurus 
maps. Overall, the observations show that, on average, the ratios of the major to minor axes 
of cores vary approximately between 1.2 to 2 with a mean value of ~ 1.6. 

Determining the exact three-dimensional shape of a core from the apparent observations 
of the-plane-of-sky is impossible, instead statistical techniques have to be applied. Several 
studies have analyzed the observations in the context of a random distribution of inclina- 
tions to infer that cores are more nearly prolate than oblate (e.g., Curry 2002, Jones and 
Basu 2002). This prolate elongation of cores may be inferred as a remnant of their origin 
in filaments (e.g., Hartmann 2002). Also, models with cores forming from turbulent flows 
predict random triaxial shapes with a slight preference for prolateness (Gammie et al. 2003, 
Li et al. 2004). Although, some works indicate a preference for prolate cores, but there 
are many studies that consistently favour oblate shapes (e.g., Jones et al. 2001, Goodwin et 
al. 2002, Tassis 2007, Offner and Krumholz 2009, Tassis et al. 2009). If strong magnetic fields 
are present then collapsing cores are expected to be oblate (e.g., Galli and Shu 1993, Basu 
and Ciolek 2004, Ciolek and Basu 2006), which could also be caused by strong rotational 
motion (e.g., Cassen and Moosman 1981, Terebey et al. 1984). 

Historically, in primary theoretical models for the collapse of cores and formation of 
star, an isothermal equation of state had been used, which as a consequence of subsonic 
communications in different parts of the cloud, an inverse square profile for density was 
appeared (Bodenheimer and Sweigart 1968). Larson (1969) and Penston (1969) were the 
first to analyze this inverse square behavior of density profile using the similarity method. 
In this case, which had been extended afterwards by Hunter (1977), one begin with a static 
cloud of constant density and follow the formation of the r~^ density profile. In the opposite 
case, Shu (1977) assumed that the density is initially in inverse square profile SIS core, and 
constructed the expansion wave collapse solution to suggest the inside-out collapse scenario. 
These two limiting solutions of Hunter and Shu may be described as fast and slow collapse. 
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respectively, and the reality may be somewhat in between (McKee and Ostriker 2007). Since 
then, a lot of asymptotic solutions and global numerical simulations have been found and 
developed in which authors have approximately considered the effects of three important 
mechanisms: turbulence, rotation and magnetic field (see, e.g., Hartmann 2009). 

In this research, we return to the basic spherical collapse problem for polytropic spheres, 
which was used as idealized SIS similarity solutions by Shu (1977). In addition, we consider 
an initial inward flows, i.e., the conditions observed in some molecular cloud cores and used 
by Fatuzzo et al. (2004) to investigate its effects on the collapse of SIS. The goal of this 
paper is to reexamine the gravitational collapse of SIS with focus on the oblateness of a 
core via the effect of rotation and magnetic field. We suggest that the centrifugal force and 
magnetic field pressure lead to oblateness of the envelope of a molecular core before reaching 
the collapsing expansion wave. Since the collapse is not in a spherical symmetric manner, 
the similarity method can not be used. Instead, we use the Adomian decomposition method 
(Adomain 1994), to solve semi-analytically the system of differential equations. For this 
purpose, the collapse of SIS using the Adomian method is given in section 2. The non- 
similar collapse of SIS is investigated in section 3 in which the oblateness of the core is also 
obtained. Finally, the section 4 devotes to summary and conclusion. 



2. Collapse of SIS by Adomian method 

The initial density of a SIS is assumed to be in the form of inverse square, and the 
collapse is assumed to be only in the radial direction. The mass continuity equation is 

where u is the radial velocity which follows the force equation, 

du du al dp , . 

dt dr p Or 

where a^ is the sound speed, and the gravitational acceleration g follows the Poisson's equa- 
tion, 

i|-(r^,) = 4.rGp. (3) 



Choosing the sound speed as the unit of velocity, and [t] = l/v^vrGp^ as the unit of time 

where pc is the density at the center of core, the basic equations (II])-(I3]) can be rewritten as 

follows 

dp dp du 2 
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du du 1 dp 

^|(^'^) = A (6) 

where ^ = r/(as[t]) is the nondimensional radius and the density and gravitational acceler- 
ation are non-dimensionalized by pc and as/[t], respectively. 

If the cloud is in hydrostatic equilibrium {u = 0), the equations (jl])-© lead to the 
Lane-Emden equation which is well-known in the theory of stellar structure (e.g., Chan- 
drasekhar 1939). The Lane-Emden equation can be solved by the Adomian decomposition 
method which is found in the appendix A. The SIS is a special case of the Lane-Emden 
equation with inverse square density which do not conform the boundary conditions. Here, 
we choose the initial density as p{^, 0) = A/^'^ where A is the overdensity parameter with 
A = 2 for hydrostatic equilibrium. The boundary condition for gravitational acceleration is 
g{0,t) = 0, and we assume that the initial velocity is inward so that u{^,0) = —u^o where 
Moo is constant. 

As mentioned in the appendix A, for using the Adomian decomposition method, equa- 
tions (jl])-© must be rewritten as 

Ltp + A^i(p,M)=0, (7) 

LtU + g + N2{p,u) = Q, (8) 

h^g - P = 0, (9) 
where Lt{o) = -^ and L^^{o) = ^^ (^^ d£ ) ^^^ operators, and 

Ni{p,u) = u-^ + p— + -pu, (10) 

^r , ^ du Idp ,_,_,, 

iV,(p.„).„_ + _^, (U) 

are the nonlinear terms of the differential equations which can be written by Adomian series 
Nj{p, u) = Y.'^=o ^n for j = 1, 2 where 



Ali^ = ^ lim 



072 

i=0 i=0 



(12) 



are the Adomian polynomials (Adomian 1994). In this way, the final solutions are given by 
the series p = J2'^=o Pn, u = Yl'^=o '^n and g = Yl'^=o 9^ which the terms can be obtained by 
the recurrence relations 

p„+i = -L;'AII\ (13) 
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Un+i = -L;\gn + A^^^), (14) 

gn+i = L^^Pn, (15) 

where L'['^{o) = j{o)dt and L'^^{o) = J ^^^ (/^^(°)t^0 ^^ ^^^ ^^^ integration operators and 
the n = terms are from initial and boundary conditions: po = p(^)O), Uq = u{C,,0) and 
go ^g{0,t). 

Using the mathematical softwares such as Maple (see, appendix B), the density and 
inward velocity are obtained as follows 

A 



P-^ 



• (A-2)t^ 2(A-2) t3 ^ (A-2)(6-2A + 3ng,) t^ , ^ f t' 



(16) 



t 



u = -Uoo - (A- 2) 



Moot , iQ-A + 2ul)e (24 - 5A + 3uy Uoo t 



1--^- + 



2 e 6 e 12 e' 



+ 



(17) 
respectively. In the case of Uoo = 0, these results reduce to the equation (19) of the well- 
known paper of Shu (1977) in which the similarity variable x is replaced by t/^ (note that 
equations (45) and (46) of Fatuzzo et al. 2004 are mistyped). According to the convergence 
problem in the series of Adomian decomposition method (as mentioned in appendix A), the 
results (IT6|) and (IT7|) are reliable in the range of t/^ < 1. Thus, the solutions flT6|) and flT7|) 
are acceptable only for the outer regions from the head of the expansion wave (i.e., C, = t), 
and the mass infall rate M = An^'^p | n | of the shells in the envelope can be determined. For 
the inner regions of expansion wave, the terms of series fITBl) and flT7|) are not convergent, 
thus, we must use the suitable methods such as piecewise-adaptive decomposition method 
(e.g., Ramos 2009) which is beyond the scope of this research. 



3. Non-similar collapse of SIS 

The cores of molecular clouds rotate and the magnetic fields are also the non-eliminating 
parts of this medium. In this section, we investigate the effect of these mechanisms on the 
dynamics of collapsing SIS core. In this case, the infall velocity of matter in the spherical 
coordinate have two components as v = f-u + 6w so that the mass continuity equation is 
expressed as 

which is non-dimensionalized according to the units given in section 2. Here, the effect of 
the centrifugal force and pressure of magnetic field, in the envelope of the core, are assumed 
to be in the normal direction of the rotational axis as follows 

F„,ag+rot = fF^ siu^ 6 + OFe sin 26, (19) 
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where for simplicity, the components are assumed to be a fraction of the local gravitational 
force as F^ = —Pd(j)/d^ and Fg = ~[5d(j)/d9 where the magnetic-rotational parameter < 
/3 < 1 indicates the importance of rotation and magnetic field, and is the gravitational 
potential which follows the Poisson's equation. 






+ 



1 d ^ . ^dct>^ 


(20) 


components as follows 




-^l-*-''-^^' 


(21) 


-^|-(-^--'||- 


(22) 



du du w du 

dw dw w dw 
'dt '^ ^~di ^ J~de '' 

Since the components of gravitational force, dcp/d^ and d(f)/d6, are zero at the center 
of the core, the boundary condition of the gravitational potential at ^ = is assumed to 
be 0(0, t) = 0. At the beginning of collapse, the components of velocity are assumed to be 
ti(r, 0) = — Uoo and w{r,0) = 0, and the initial density is assumed as the density of SIS: 
p(^, 0) = A/^^ We rewrite the the equations ([H]), ([20D, ([2l]) and ([22]) in the Adomian form 
as follows 

Ltp + Ni{p,u,w) = 0, 



^C50 + 



1 92 



^2 ^^2 

Ltu+ (1 -/3sin2 



+ 



cot 6'( 



P 



d^ 



^2 de 

+ N2{p,u,w 



0, 



Ltw + (1 - /3sm29)-^ + N:i{p,u,w) = 0, 

respectively, where the nonlinear terms are 

, ^ , , dp du 2 w dp p dw cot 6 

Ndp,n,u.) ^u^^+p^ + -pu + -_ + -_ + _ 

^^ , , du w du 1 dp 

iV,(p,n,^)^.- + -- + --, 

dw w dw 1 dp 
Nsip,u,w)^u- + -- + --. 

In this way, by appointment the Adomian polynomials, 

1 



pw, 



Ay^ 



lim 

n! A^o 



d\^ ' 



^PiX\^UiX\^w,X' 

J=0 i=0 i=0 



(23) 
(24) 

(25) 

(26) 

(27) 

(28) 
(29) 

(30) 
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for j = 1, 2, 3, we receive to the recurrence relations as follows 



■>n+l 



Un+1 



Wn+1 



-L^i 



'L;' 



1 92 



+ 



cot 6' I 



^2 QQ2 ^2 QQ 



Pr. 



/3sin2 



d^ 



^ + Al^^ 



-L 



-1 



(1-/3 sin 2e)i^ + Af 



^ 90 



(31) 
(32) 

(33) 

(34) 



where the n = terms are given by the initial and boundary conditions. Thus, we can obtain 
the terms of the series p = Y.n=oPn, u = J2n=o^n, w = Y.n=o^n and = Y.n=o4>n with a 
mathematical software such as Maple (see, appendix B). Here, we turn our attention to the 
mass infall rate, M = An^'^p | m | , at the envelope of the core (^ > t) where the Adomian 
decomposition method is reliable. The result is 



M = AttAuoo + 47rA(A - 2 - A/3 sin^ 9)-- AttA{A - 2 - A/3 sin^ 1 



t' 



]u. 



°°t2 +^ V f3 



t' 



(35) 



which its values for A = 2.1, Uao = 0.1 and /3 = 0.1, at time t = 1/6 are shown in FigJT] 

Since the infall rates at the polar and equatorial regions are different, an initial spherical 
shell at the outer region of expansion wave, will be spheroid by time, as schematically shown 
in Fig. HJ In the first order approximation of mass infall rate (i.e., in the order of t/C,), ratio 
of major to minor axes of the spheroid can be approximated as 



b 
a 



^ - u^t - |M(e=^/2)t^ 
e - Moot - lM(e=o)t^ 



(36) 



Inserting M ~ 47rA(A — 2 — A/3 sin^ ^)/^, which is obtained from equation ( l35l) . we have 



b_ 
a 



l-Moo7-27rA(A-2-A/3)|^ 



l-Moof-27rA(A-2)| 
which is depicted in Fig. |3]at time t = 1/6. 



(37) 



4. Summary and conclusions 

Stars have been formed from collapse of cores in the molecular clouds. A basic standard 
model for collapse of cores assumes a SIS, in which the collapse occurs inside-out accompanied 



with an expansion wave. We know that not only the molecular cores rotate, but also the 
magnetic fields affect on their dynamics. If we assume that the effect of centrifugal force 
and magnetic pressure are in the normal direction of the rotational axis, we expect that the 
shape of the envelope of core (outer regions from the head of expansion wave) be modified 
to a spheroid. In this case which the collapse is in a non-symmetric manner, we used the 
Adomian decomposition method to solve the differential equations. In the appendix A, we 
solved the well-known Lane-Emden equation to find that the Adomian method is applicable. 
Then, in section 2, the collapse of a SIS core is investigated. The results show that the 
Adomian decomposition method presents convergent solutions only for the regions outer 
than the head of the expansion wave (i.e., envelope). 

In section 3, non-similar collapse of a SIS core, which is affected by the centrifugal force 
and magnetic pressure, is solved by the Adomian decomposition method. The mass infall 
rate of the shells in the envelope of the core is obtained and depicted in FigJT] with typical 
values of the parameters. According to Fig. [TJ^a), the mass infall rate in the polar regions 
(near the rotational axis) decreases with increasing of radius, while in equatorial regions, the 
mass infall rate in the regions far from the head of expansion wave is greater than regions 
near to it. Fig. Wljo) shows that the mass infall rate in equatorial regions are less than 
polar regions, thus, before reaching the expansion wave, shape of the envelope converts to a 
spheroid as shown schematically in the Fig. [2l 

In this way, according to mass infall rate, we determined the ratio of major to minor axes 
of the spheroid, and its values are shown in the Fig. Elwith typical values of the parameters. 
We see that the ratio of major to minor axes of spheroid do not strongly depend on the initial 
inward velocity Mqo, while its value depends on the overdensity parameter A and magnetic- 
rotational parameter /3. Thus, the results show that the magnetic field and rotation lead to 
non-symmetric collapse so that before reaching the head of expansion wave, the apparent 
shape of the envelope will be converted to the triaxial oblate spheroids. Also, the results 
show that the ratio of major to minor axes of the spheroids can reach to the values which 
are consistent with the observed maps of the cores in the molecular clouds such as Taurus 
and Perseus. 
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A. Solving the Lane-Emden equation by Adomian decomposition method 



A spherical polytropic gas, in hydrostatic equilibrium, satisfies the equation 

ldp_ GMjr) 
pdr ^2 ' 



(Al) 



where M{r) is the mass inside the radius r and the polytropic pressure is p = /tp'", which 
in the isothermal case, we have n = al and F = 1, where a^ is the isothermal sound speed. 
Substituting p = pcexp{—ip) into equation (lAip . we have 






2^(i-r)^ 



Ce 



dip 



(A2) 



1/2 



where ^ = ( ^""/?2 ) is the nondimensional radius. In the isothermal case, the equation 
( ]A2p reduces to the well-known Lane-Emden equation. The boundary conditions are 4'(^=o) = 
and {dip/d^)f^Q = 0. Thus, the equation ( JA2I) can straightforwardly be solved by numerical 
methods such as fourth-order Runge-Kutta. 



where Lc 



Here, we solve the equation (IA2I) by the Adomian decomposition method, then we 
compare the result with the numerical method. For this purpose, we rewrite the equation 
l21) as follows 

L^^^ + A(^) = 0, (A3) 

) = ^2^ (^^ df ) ^^ ^^^ operator, and nonlinear terms of the differential equation 

are assembled in the function N{iIj) = (1 — F) ( -^ j — e^^~'^'>'^. The basis of the Adomian 
decomposition method is to replace the function ■?/' by a series ip = S^o"^"' ^^^ ^^^ ^°^" 
linear term N{il)) by a Taylor expansion series N{iIj) = ^.^g An where An are Adomian 
polynomials 

Ao = A(^o), 



V'=V'o 



-^-'l)..4^? 



d^^ 

He 



(A4) 



■>P=i'o 



and so on (Adomian 1994), which can be expressed in the general form 



An = —r lim 

nl x^o 



dX^ 



-nIj^^pa^ 



j=o 



(A5) 
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where A is only a dummy variable which is inserted to recover the equations flA4|) . 
In this way, the equation (IA2p leads to a recurrence relation 

V^„+i = -L7M,, (A6) 



where L^l{o) = f (.^^ {f ^^{°)d^) d^ is the integration operator and ipo = 4'(^=o)- Thus, 
we can determine the ipnS with help of any mathematical softwares such as Maple (see, 
appendix B). The result is a series as 



^ = 02^2 + 04^^ + 06^^ + ... 



(A7) 



where the coefficient are 



«2 = +7:7 



a. = --i 



ae = +: 



4(2-r)2 + 5(2-r)(l 



r) + 3(i 



r)^ 



1 

7!^ 



(A8) 



and so on. In isothermal case (F = 1), this result reduces to the result of Wazwaz (2001). 

e~^, in the isothermal case, which 



In Fig. m the ratio of density to central density, — 



is obtained from three terms of the series (lASp is compared to the numerical results, which 
are obtained from the fourth-order Runge-Kutta method. As can be seen, the Adomian 
decomposition method gives reliable solutions for ^ < 2. This is because of using the 
Taylor expansion which destroys the convergence of the series at large ^ (Liao 2003, Singh 
et al. 2009). Substituting the appropriate values of the typical molecular clouds, we have 



e 



7pc) 



0.03pc(^)^^^(,P^) 



-1/2' 



(A9) 



where n^ is the number density at the center of the core. Thus, in the hydrostatic equi- 
librium of a typical molecular cloud core, the results which are obtained by the Adomian 
decomposition method, are usable only for radii less than ~ 0.06pc. 



B. Maple programs 

The Maple program for the Lane-Emden equation is as follows: 

> psi[0] := 0: 

> n:= 0: 

> N:= proc (psi) options operator, arrow; 



-li- 
d-Gamma) * (dif f (psi , xi) ) ~2-exp ( (Gamma-2) *psi) 
end proc : 

> psi[l]:= -int(xi"(-2)*int( xi~2*N(psi [0] ) ,xi) ,xi) : 

> for n from 1 to 7 do 

N:= proc (psi) options operator, arrow; 

(1-Gamma) * (dif f (psi ,xi) ) "2-exp( (GcLmma-2) *psi) 
end proc: 
psisum:= sum(psi [m] *lambda[n] "m,m=0. .n): 
A[n] := (diff (N(psisum) ,$(lainbda[n] ,n)))/f actorial(n) 
lambda [n] := 0: 

psi[n+l] := -int(xi"(-2)*int(xi"2*A[n] ,xi) ,xi) : 
end do : 

> psitotal:= collect (sum(psi [m] ,m=0. .n) ,xi) ; 

The Maple program for the collapse of SIS is as follows: 

> g[0]:= 0: 

> rho[0]:= Lambda/xi"2: 

> u[0] := -u [infinity] : 

> n:= 0: 

> Nl:= proc (rho,u) options operator, arrow; 

u*dif f (rho , xi) +rho*dif f (u , xi) +2*u*rho/xi 
end proc : 

> N2:= proc (rho,u) options operator, arrow; 

u*diff (u,xi)+rho"(-l)*diff (rho,xi) 
end proc : 

> g[l]:= xi"(-2)*int(xi~2*rho[0] ,xi) : 

> rho[l]:= -int(Nl(rho[0] ,u[0]),t) : 

> u[l]:= -int(g[0]+N2(rho[0] ,u[0]),t): 

> for n from 1 to 7 do 

Nl:= proc (rho,u) options operator, arrow; 

u*dif f (rho , xi) +rho*dif f (u , xi) +2*u*rho/xi 
end proc: 
N2:= proc (rho,u) options operator, arrow; 
u*diff (u,xi)+rho~(-l)*diff (rho,xi) 
end proc : 
rhosum:= sum(rho[m] *lambda[n] "m,m=0. .n): 
usum:= sum(u[m] *lambda[n] "m,m=0. .n) : 



-12- 



Al[n] := (diff (Nl(rhosum,usum) ,$(lambda[n] ,n)))/f actorial(n) : 
A2[n] := (dif f (N2(rhosum,usum) ,$(lambda[n] ,n)))/f actorial(n) : 
lambda [n] := 0: 

g[n+l] := xi"(-2)*int(xi"2*rho[n] ,xi) : 
rho[n+l]:= -int(Al [n] ,t) : 
u[n+l]:= -int(g[n]+A2[n] ,t) : 
end do : 

> rhototal:= collect (sum(rho [m] ,m=0. .n) ,t) ; 

> utotal:= collect(sum(u[m] ,m=0. .n) ,t) ; 

> gtotal:= collect(sum(g[m] ,m=0. .n) ,t) ; 

The Maple program for the non-similar collapse of SIS is as follows: 

> Phi [0] : = : 

> rho[0]:= Lambda/xi~2: 

> u[0] := -u [infinity] : 

> w [0] : = : 

> n:= 0: 

> Nl:= proc (rho,u,w) options operator, arrow; 

u*dif f (rho , xi) +rho*dif f (u , xi) +2*u*rho/xi+rho/xi*dif f (w , theta) 
+w/xi*dif f (rho , theta) +cot (theta) *rho*w/xi 
end proc : 

> N2:= proc (rho,u,w) options operator, arrow; 

u*diff (u,xi)+w/xi*diff (u,theta)+rho"(-l)*diff (rho,xi) 
end proc : 

> N3:= proc (rho,u,w) options operator, arrow; 

u*dif f (w , xi) +w/xi*dif f (w , theta) +rho" (-1) *dif f (rho , theta) /xi 
end proc: 

> Phi[l]:= -int(xi"(-2)*int(xi~2*(diff (Phi[0] , theta, theta)/xi"2 

+cot (theta) *diff (Phi [0] , theta) /xi~2-rho [0] ) ,xi) ,xi) : 

> rho[l]:= -int(Nl(rho[0] ,u[0] ,w[0]),t): 

> u[l]:= -int((l-beta*(sin(theta))~2)*diff (Phi[0] ,xi) 

+N2(rho[0],u[0] ,w[0]),t): 

> w[l]:= -int((l-beta*sin(2*theta))*diff (Phi [0] , theta) /xi 

+N3(rho[0] ,u[0] ,w[0]),t) : 

> for n from 1 to 7 do 

Nl:= proc (rho,u,w) options operator, arrow; 

u*dif f (rho , xi) +rho*dif f (u , xi) +2*u*rho/xi+rho/xi*dif f (w , theta) 
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+w/xi*dif f (rho , theta) +cot (theta) *rho*w/xi 
end proc : 
N2:= proc (rho,u,w) options operator, arrow; 

u*dif f (u , xi) +w/xi*dif f (u , theta) +rho~ (- 1) *dif f (rho , xi) 
end proc : 
N3:= proc (rho,u,w) options operator, arrow; 

u*dif f (w , xi) +w/xi*dif f (w , theta) +rho~ (- 1) *dif f (rho , theta) /xi 
end proc : 
rhosum:=sum(rho [m] *lambda[n] "m,m=0. .n) : 
usum:=sum(u[m] *lambda[n] ~m,m=0. .n) : 
wsum:=sum(w[ni] *lambda[n] ~m,m=0. .n) : 

Al [n] := (dif f (Nl(rhosuni,usum,wsuni) ,$ (lambda [n] , n) )) /factorial (n) 
A2[n] := (dif f (N2(rhosuni,usum,wsuni) ,$ (lambda [n] ,n)))/f actorial(n) 
A3[n] := (dif f (N3(rhosuni,usum,wsuni) ,$ (lambda [n] ,n)))/f actorial(n) 
A4[n] := (dif f(N4(rhosum, usum,wsum) ,$ (lambda [n] ,n) )) /factorial (n) 
lambda [n] :=0: 

Phi [n+1] : =-int (xi" (-2) *int (xi"2* (dif f (Phi [n] , theta, theta) /xi~2 
+cot (theta) *diff (Phi [n] , theta) /xi ~2-rho [n] ) ,xi) ,xi) : 
rho [n+1] :=-int(Al[n] ,t) : 

u[n+l] : =-int((l-beta* (sin (theta) )~2)*diff (Phi [n] ,xi)+A2[n] ,t) : 
w[n+l] :=-int((l-beta*sin(2*theta))*diff (Phi[n] , theta) /xi+A3 [n] ,t) 
end do : 

> rhototal:=collect(sum(rho [m] ,m=0. .n) ,t) : 

> utotal:=collect(sum(u[m] ,m=0. .n) ,t) : 

> wtotal:=collect(sum(w[m] ,m=0. .n) ,t) : 

> phitotal:=collect(sum(Phi [m] ,m=0. .n) ,t) : 

> mdot:=collect(4*Pi*xi~2*(rhototal)*(-utotal) ,t) ; 
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Fig. 1. — Mass infall rate versus (a) nondimensional length ^, and (b) polar angle 6, with 
A = 2.1, Moo = 0.1 and (3 = 0.1, at time t = 1/6. 
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Fig. 2. — Schematic diagram for mass infall rate which leads to the formation of a spheroidal 
shell from an initial spherical shell. 
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Fig. 3. — The ratio of major to minor axes of spheroidal shells of a core, at time t = 1/6, 
versus (a) magnetic-rotational parameter /3 with A = 2.1 and Uoo = 0.1, (b) overdensity 
parameter A with /3 = 0.1 and Uoo = 0.1, and (c) inward initial velocity Uoo with A = 2.1 
and (3 = 0.1. 
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Fig. 4. — Ratio of density to the central density of an isothermal spherical core in the 
hydrostatic equilibrium. 



